#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <assert.h>

/*******************************************************************************
 *
 * utilities
 *
 ******************************************************************************/

double alea()
{
    static const double fac = 1.0 / ( ((double)RAND_MAX) + 1.0 );
    return fac * ((double)rand()+0.5);
}


void InitializeSimulation()
{
    int i=0;
    srand( time(NULL) );
    for(i=0;i<1024;++i)
        (void) rand();
}

double hmax             = 1e-4;

double rk4( double (*f)(double), double t1, double t2 )
{
    int    i;
    double s     = 0;
    double half  = 0;
    double h     = hmax;
    double h6    = 0;
    int    nstep = ceil(fabs(t2-t1)/fabs(h));
    if(nstep<1) nstep=1;
    h    = (t2-t1)/nstep;
    half = h*0.5;
    h6   = h/6.0;
    
    for(i=0;i<nstep;++i)
    {
        const double t  = t1 + i*h;
        const double k1 = f(t);
        const double k2 = f(t+half);
        const double k3 = f(t+half);
        const double k4 = f(t+h);
        s += h6*(k1+(k2+k2)+(k3+k3)+k4);
    }
	
    
    return s;
}



/*******************************************************************************
 *
 * constants
 *
 ******************************************************************************/

double propagation_time = 0.1;
double deltaP           = 0.1;
double omega            = 1.0;
double beta             = 0.1;
double zeta             = 1.0;
double P_star           =-10.0;

/*******************************************************************************
 *
 * Description of one cavity
 *
 ******************************************************************************/

#define IS_LIQUID 0
#define IS_BUBBLE 1

#define MAX_NEIGHBORS 2

typedef struct Cavity_
{
    int    status; /*!< what kind of cavity */
    double t_cav;  /*!< cavitation time     */
    int    pos;    /*!< the x index         */
    double lam;    /*!< the intensity       */
    
    int             neighbors;               /*!< 0..MAX_NEIGHBORS      */
    struct Cavity_ *neighbor[MAX_NEIGHBORS]; /*!< pointers to neighbors */
    int             tprograg[MAX_NEIGHBORS]; /*!< propagation times     */
    
}  Cavity;

int     nx      = 1000;
Cavity *cavity  = NULL;
int     bubbles = 0;

/*******************************************************************************
 *
 * Memory
 *
 ******************************************************************************/
void AllocateCavities()
{
    cavity = (Cavity *) calloc( nx, sizeof(Cavity) );
    if(!cavity)
    {
        perror("CavityAllocate");
        exit(1);
    }
}

void DeallocateCavities()
{
    if(cavity)
    {
        free(cavity);
        cavity = NULL;
    }
}

/*******************************************************************************
 *
 * Initialize
 *
 ******************************************************************************/
void InitializeCavities()
{
    int i;
    for(i=0;i<nx;++i)
    {
        Cavity *c = &cavity[i];
        /* erase everything */
        memset(c,0,sizeof(Cavity));
        
        /* set up things */
        c->status = IS_LIQUID;
        c->pos    = i;
        c->t_cav  = -1;
    }
    bubbles = 0;
}

/*******************************************************************************
 *
 * modify one cavity
 *
 ******************************************************************************/
void CavityToBubble(Cavity *c,double t)
{
    assert(NULL!=c);
    assert(t>=0);
    assert(bubbles<nx);
    c->status = IS_BUBBLE;
    c->t_cav  = t;
    ++bubbles;
}



/*******************************************************************************
 *
 * detect neigbhors
 *
 ******************************************************************************/
void MakeNeighborOf( Cavity *c, Cavity *p )
{
    assert(c->neighbors<MAX_NEIGHBORS);
    assert(IS_LIQUID==c->status);
    assert(IS_BUBBLE==p->status);
    c->neighbor[c->neighbors] = p;
    c->tprograg[c->neighbors] = propagation_time * abs(c->pos-p->pos);
    c->neighbors++;
}

void DetectNeighborsOf(Cavity *c)
{
    int i;
    assert(NULL!=c);
    assert(IS_LIQUID==c->status);
    c->neighbors = 0;
    
    /* check left */
    for(i=c->pos-1;i>=0;--i)
    {
        Cavity *p = &cavity[i];
        if(IS_BUBBLE==p->status)
        {
            MakeNeighborOf(c,p);
            break;
        }
    }
    
    /* check right */
    for(i=c->pos+1;i<nx;++i)
    {
        Cavity *p = &cavity[i];
        if(IS_BUBBLE==p->status)
        {
            MakeNeighborOf(c,p);
            break;
        }
    }
    
}

void DetectAllNeighbors()
{
    
    int i;
    for(i=0;i<nx;++i)
    {
        Cavity *c = &cavity[i];
        if(IS_LIQUID==c->status)
        {
            DetectNeighborsOf(c);
        }
    }
}

/*******************************************************************************
 *
 * Physics: compute one BUBBLE cavity superpressure
 *
 ******************************************************************************/
double SuperPressureOf(const Cavity *c, double t)
{
    assert(NULL!=c);
    assert(IS_BUBBLE==c->status);
    if(t<=c->t_cav)
    {
        return 0;
    }
    else
    {
        const double tt = t-c->t_cav;
        /* TODO: put the physics here!!! */
        //return 0.1*sin(tt);
        return deltaP*(1.0-cos(omega*tt)*exp(-beta*tt));
    }
}

/*******************************************************************************
 *
 * Physics: compute one LIQUID cavity total superpressure
 *
 ******************************************************************************/
double ComputeIntensityOf(Cavity *c, double t)
{
    /* at first: no superpressure */
    double DeltaP = 0;
    int    i;
    assert(IS_LIQUID==c->status);
    
    /* collect superpressure of the neighbos */
    for(i=0;i<c->neighbors;++i)
    {
        const Cavity *p     = c->neighbor[i];
        const double  delay = c->tprograg[i];
        DeltaP += SuperPressureOf(p,t-delay);
    }
    
    /* TODO: convert to intensity */
    c->lam = exp(-zeta/((P_star+DeltaP)*(P_star+DeltaP)));
    
    return c->lam;
}

/*******************************************************************************
 *
 * Physics: compute sum of all intensities, function to integrate
 * WARNING: Neighbors MUST be detected
 *
 ******************************************************************************/
double ComputeIntensity(double t)
{
    double ans = 0;
    int    i;
    for(i=0;i<nx;++i)
    {
        Cavity *c = &cavity[i];
        if(IS_LIQUID==c->status)
        {
            ans += ComputeIntensityOf(c,t);
        }
    }
    return ans;
}

/*******************************************************************************
 *
 * Find next time
 * WARNING: Neighbors MUST be detected
 *
 ******************************************************************************/
double FindNextTime(double t)
{
    int   iter     = 0;
    const double U = -log( alea() );
    assert(bubbles<nx);
    
    /* Initialize */
    double t_last = t;
    double lambda = ComputeIntensity(t_last);
    double Lambda = rk4(ComputeIntensity,t_last,t=U/lambda);
	
    
    /* Newton Algorithm */
    /* TODO: write a convergence test */
    while(1)
    {
        
        const double residue = U-Lambda;
        const double dt      = residue/(lambda=ComputeIntensity(t));
        t_last = t;
        t+=dt;
        ++iter;
        /* fprintf(stderr,"residue=%g,t=%g\n", residue,t); */
        if(iter>=10)
            break;
        
        /* Increase Lambda from t -> t+dt */
        Lambda += rk4(ComputeIntensity,t_last,t);
        
    }
    
    return t;
}


/*******************************************************************************
 *
 * make one random bubble at cavitation time t
 *
 ******************************************************************************/
void MakeBubble(double t)
{
    const double Lambda = ComputeIntensity(t);
    const double Level  = Lambda * alea();
    double       sum    = 0;
    int          i;
    assert(bubbles<nx);
    
    for(i=0;i<nx;++i)
    {
        Cavity *c = &cavity[i];
        if(IS_LIQUID==c->status)
        {
            sum += c->lam;
            if( Level <= sum )
            {
                /* found ! */
                CavityToBubble(c,t);
                return;
            }
        }
    }
    
    /* should NEVER get here */
    fprintf(stderr,"failed to find one cavity for cavitation!\n");
    exit(2);
}



/*******************************************************************************
 *
 * debugging
 *
 ******************************************************************************/

void Display()
{
    int i;
    for(i=0;i<nx;++i)
    {
        switch(cavity[i].status)
        {
            case IS_LIQUID:
                fprintf(stderr,"#");
                break;
                
            case IS_BUBBLE:
                fprintf(stderr,"O");
                break;
                
            default:;
        }
    }
    fprintf(stderr,"\r");
    fflush(stderr);
}

//______________________________________________________________________________
//
//
// Main
//
//______________________________________________________________________________

int main(int argc, char *argv[])
{
    double t =0;
    FILE   *fp = NULL;
    InitializeSimulation();
    AllocateCavities();
    
    InitializeCavities();
    t=0;
    fp = fopen("output.dat","wb");
    fprintf(fp,"%g 1\n",t);
    Display();
    while(bubbles<nx)
    {
        DetectAllNeighbors();
        t = FindNextTime(t);
        MakeBubble(t);
        Display();
        printf("\n");
        fprintf(fp,"%g %g\n",t,(nx-bubbles)/(double)nx);
    }
    fprintf(stderr,"\n");
    fclose(fp);
    DeallocateCavities();
    return 0;
}

